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Abstract 

The fluctuation-dissipation (F-D) theorem is a fundamental result for 
systems near thermodynamic equilibrium, and justifies studies between 
microscopic and macroscopic properties. It states that the nonequilibrium 
relaxation dynamics is related to the spontaneous fiuctuation at equilib- 
rium. Most processes in Nature are out of equilibrium, for which we have 
limited theory. Common wisdom believes the F-D theorem is violated in 
general for systems far from equilibrium. Recently we show that dynam- 
ics of a dissipative system described by stochastic differential equations 
can be mapped to that of a thermostated Hamiltonian system, with a 
nonequilibrium steady state of the former corresponding to the equilib- 
rium state of the latter. Her we derived the corresponding F-D theorem, 
and tested with several examples. We suggest further studies exploiting 
the analogy between a general dissipative system appearing in various 
science branches and a Hamiltonian system. Especially we discussed the 
implications of this work on biological network studies. 



It is ubiquitous to observe a system at a state invariant with time (with the 
approximation that the relevant constraining parameters changes much slower 
than the time scale under interest). It can be a thermodynamic equilibrium 
state, or more likely a nonequilibrium steady-state. Some examples are homeo- 
static states of living organisms, a stable eco- or financial system. Quite often it 
is important to know how a system initially subject to a perturbation relaxes to 
a steady state (including the equilibrium state) after removal of the perturba- 
tion. The fluctuation-dissipation theorem states that the relaxation dynamics 
for a process close to equilibrium is related to the spontaneous fluctuation at 
equilibrium. Originally formulated by Nyquist in 1928 [T], and first proved by 
Callen and Welton in 1951 [5], the F-D theorem is related to many important 
results in statistical physics. Examples are the Einstein- Smoluchowski relation 
on the diffusion constant and drag coefficient [3j|4], Onsager's regression hy- 
pothesis |BJ H] , and the linear response theory [7] . The F-D relation also has 
practical importance. It allows deducing nonequilibrium dynamics from equi- 
librium measurements, and justifies the relation between macroscopic dynamics 
and microscopic level simulations, e.g., calculating the diffusion constant. In 
recent years, fiuctuation theories of nonquilibrium processes, especially for sys- 
tems far from equilibrium received great attention [51 [HI UHl [HI HI] ■ 

On studying problems from physics, chemistry, cellular biology, ecology, en- 
gineering, finance, and many other fields, the following form of equations are 
widely used [Il[Til[IS], 

m 

dxi/dt ^ Gi(x) + ^5ij(x)Cj(t), i = 1, • • • ,n. (1) 

In general m and n may be different, Cj(^) ^-re temporally uncorrelated, statisti- 
cally independent Gaussian white noise with the averages satisfying < Q {t)Cj (t) > — 
Sjji5{t — t), g(x) is related to the nxn diffusion matrix gg^ = 2D//3, where the 
superscript T refers to transpose. For a physical system /3 is the inverse of the 
Boltzmann's constant multiplying temperature. For a non-physical system, one 
can define an effective temperature relating to /3. Recently we proved that there 
exists a mapping between a stochastic dissipative system described by Eqn. [1] 
and a thermostated Hamiltonian system [T7]. The mapping allows many re- 
sults from equilibrium statistical physics directly applicable to nonequilibrium 
processes. Specifically in this work, we will derive the F-D relation applica- 
ble to processes far from equilibrium. An F-D relation, if exists, would allow 
predicting the relaxation dynamics to a steady-state through measurements of 
steady-state fluctuations. The latter are in general easier to measure. 



1 



1 Theory 



1.1 Existence of mapping between linear stochastic dissi- 
pative systems and Hamiltonian systems 

First we briefly summarize the main result of [17 . The mapping is based a 
seminal work of Ao, which shows that one can always construct a symmetric 
matrix S and an anti-symmetric one T, and transform Eqn. [T]into, |18) 



(S + T)(G(x)+g(x)C(t)) 
-Vx0(x)+g'(x)C(t) 



(2) 



where (/> is a scalar function corresponding to the potential function in a Hamil- 
tonian system satisfying {d x d(j))ij = {didj — djdi)<j) — 0, and g'g'^ = 2S//3. 
Then S and T are uniquely determined by 



d X [(MG(x)] 0, (M)-i + (M)- 



2gg' 



(3) 



where M = S -f- T, with proper choice of the boundary conditions. In [17], we 
first demonstrated that starting with the corresponding Fokker-Planck equa- 
tions, one can derive the transformation matrix more transparently following a 
standard procedure used by Graham and by Eyink et al. previously [lH [20] . 
Then we showed that one can map the dynamics described by Eqn. [2] to a 
Hamiltonian system in the zero mass limit. 



A(x))^ 



2m 



No. 

■E 

a = l 



N 

E 



1 2 



H 



(g,, -a,(x)/(ViVc^^,)) 



(4) 



dxi 



The last 



The term A is a vector potential satisfying Tij = 

term is a type of bath Hamiltonian discussed by Zwanzig [3T]. The Hamilto- 
nian mathematically corresponds to Dirac's constrained Hamiltonian [21]. It 
describes a massless particle, coupled to a set of harmonic oscillators, mov- 
ing in a hypothetical n-dimensional conservative scalar potential and magnetic 
(the vector potential) field. The steady state distribution is thus given by the 
Boltzmann distribution of the Hamiltonian system, 



Pss{x) 



exp( 



/dp<iYexp(-/3i/) 
J dxdpdY exp(— J dxexp(- 



(5) 



where Y represents all the bath variables. Eqn. [5] is also conjectured by Ao 
[T5] , and a general proof is given in [53] . 
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1.2 General fluctuation-dissipation theorem 

Let's denote the steady-state ensemble average of a generic dynamic quantity 
O as, 

< O >= / dxOp,,(x) (6) 



Consider a system initially at a steady state defined by Eqn. [T] with an ex- 
tra infinitesimal perturbation ^G(A), where A refers to the parameters be- 
ing perturbed. The corresponding perturbation Hamiltonian term is given by 
V{6H) = -(M'G' - MG) « -MSG - SMG, where SM is the variation of M 
due to the perturbation. At time SG is removed, and the system relaxes to 
the steady-state of G. Or the Hamiltonian of the corresponding mapped system 
changes from H' = H + SH = H + ^ ■ SX. Then the nonequilibrium relaxation 
dynamics of O follows 

- , , _ f dxdpO(t)exp(-/3H') 

0{t)- < O r/^ ^ /^ff/A < O > 

J c(xdpexp(— pii') 



-/?[<o(t)^^)>-<o><|>].a(o) 



(7) 



In the above derivation, we used the stationary property of equilibrium ensemble 
average, so < 0{t) >—< 0(0) >=< O >. This is the generalized F-D relation 
for systems obeying or violating detailed balance, which states that nonequilib- 
rium relaxation dynamics can be predicted from steady-state fluctuations. If 
one relates the relaxation function to the linear response function, 

0^{ty<0^>^f dt'x„{t-t')5\j{t') (8) 

"'to 

one obtains the differential form of the F-D relations, 

(9) 



dXj 



Eqns [7] and [H are the central results of this work. These results are actually 
mathematically trivial with the replacement cj) — —\n pss [H] • The mapping, 
however, provides direct connection between pss and the Langevin equations 
(see especially case 4 in the next section), and allows unified treatment for 
systems obeying or violating detailed balance. One can generalize the results 
discussed in this work to higher orders oi SG. 

One type of perturbation of special interest is a system coordinate cou- 
ples to some constant external force linearly, SH = f • x, which corresponds to 
SG = — M~-'-f with — 0. Notice that SG is in general a nonlinear func- 
tion of x. This situation has been previously discussed by Graham, and by 
Eyink et at [191 HO]. Under this special type of perturbation, all the famil- 
iar results obtained on studying relaxations near an equilibrium state follow 
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pi] . One can define a response function x{t, t') so x{t)~ < x >= dt'x{t — 
i')f(t') + 0(||/|p). The function x{t—t') is stationary and satisfies the Kramers- 
Kronig relations. If f is time varying with a monochromic frequency, f = 
Re[fi^ exp{iijjt) + f* exp(— iwt)], the system absorbs "energy", with the absorp- 
tion spectrum abs{uj) cx dti^ < (xf — < x > (xq— < x >)'^ > f* cos{ujt) 

m- 



2 Special cases and numerical tests 

Here we wih show that several versions of the previously derived generalized 
F-D relations are special cases of the above results. We will also consider sev- 
eral pedagogical examples, especially biochemical networks, to demonstrate the 
validity of our results. 

Case 1: The first one is an analytically solvable irreversible linear chemical 
network (Fig. la). The system is initially at a state with kg + Sk. Then at time 
it changes to ko, i.e., the infiow flux varies. The system relaxes to a new steady- 
state. The perturbation is 5Q = ((5fc, 0,0), and the corresponding Hamiltonian 
perturbation term is 5H — — (a;i, 2:2, X3)M(5G. The dynamic equations are. 




(10) 



For simplicity of the following discussions, we write the above equation in the 
form dx/dt = Kx + h + z((t). The predicted relaxation function through the 
fluctuation-dissipation relation is given by, 

{x^{t)~ <x,>)fd oiC-M-SG (11) 

where {Cij{t)) —< {xi{t)— < Xi >)[xj{Q)— < xj >) >. We also tested an 
incomplete relaxation function by setting M as an identity matrix, 

(x^ (i) <^ Xi ^^incomplete C • SG (12) 

Analytical solutions are given in the appendix. Fig. 1 shows that Axi (t) calcu- 
lated from the F-D relation (Eqn [Tl]) reproduces the exact result given by Eqn 
[2T] (solid line). However, without the term M, the incomplete F-D relation Eqn 
1121 predicts an initial increase of X2 (and similar for 0:3). Physically, sponta- 
neous fluctuation of xi is anti-correlated with those of X2 and X3 concentrations 
through reactions 1 and 3. The incomplete F-D relation erroneously attributes 
this mechanism to the relaxation dynamics. 

In the second example, we still use the network shown in Fig. la, but assume 
that all the reactions except reaction follow the Michaelis-Menten kinetics, e.g., 

dxi viXi V3X1 



dt Ki + xi K3 + xi 



zCi{t) (13) 
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and similar expressions for xi and xa.The relaxation functions are simulated 
numerically. Appendix B gives computational details. First the perturbation is 
again on k^. Fig. 2 shows that the relaxation dynamics at small perturbations 
can be well reproduced by the linear FD theory. Large deviation at larger 
perturbation amplitudes may be due to breakdown of the linear approximation 
for either G or in Eqn. [T] Alternatively, we added perturbation on v\. When v\ 
decreases from an initial larger value, x\ increases due to decreased consumption 
through reaction 1; X2 first decreases, then recovers to a new steady state value 
after the flux from xj, flows in ; X3 increases since the flux of the competing 
x\ — >■ xi pathway decreases, while the incoming flux to X\ is unchanged. Fig. 3 
shows that the FD theory reproduces all these behaviors. We found deviation 
only at rather large perturbation values. This example demonstrates that the 
linear approximation works for more realistic models. 

case 2: If one takes = If ., Eqn. E] gives the relation derived by Prost et 
at from the Hatano-Sasa equality recently [13] . 

case 3: Consider a particle moving along a one-dimensional periodic po- 
tential V{x + L) = V{x) and under a constant force F. The corresponding 
Langevin equation is 

dx/dt^ d^Veff + V2l5c{t) (14) 

with 14// = V ~\- Fx. The "steady state" distribution (projected to and nor- 
malized within a range of L) is , 

Pss{x)--N f dye[(-f^"("'"(^+^)-'''"^)-^^)/^l (15) 
Jo 

where M is the normalization factor. Then Eqn. [7] can be rewritten as. 



Oit)- < O ^ 



OA 



<0>< — ^if^ > 
oX 



- j^[< OB > - < O >< B >] (16) 



where B = dy ^^-ii^^^+y'^ e^iHo (sin(x+j;)-sin x)-Fv)/d] ^ rpj^g second term in the 
right hand side of the above equation vanishes when F — 0, and so the relation 
reduces to the F-D relation near equilibrium. In the special case O and ^^"gj^^^^ 
are conjugate variables, and the motion is along a circular track, it remains to 
show that the current result is equivalent to what derived by Chetrite et al. 
[27], and experimentally tested by Solano et al [28j . 

case 4: To further test the F-D model in nonlinear cases, we considered a 
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model studied by Zhu et al. 



m 



d_ 

di 



( 



) 



(S + T)-iVV' + C 



(17) 



T 



S 



jq'-l? ( 1 

(q2 _ 1)2 + 1 1^ 

(g - 1 + e)q^ / 

(g2 - 1)2 + 1 1 




< CO > 



(18) 



In the last expression, the diffusion matrix D is given by the symmetric part 
of (S + T)~^, and we used /3 — 20. Compared to what used by Zhu et ai, we 
added an extra term e = 0.01 to avoid singularity at g = 1 in the introduced 
perturbation, SG = (S + T)~^f, even with elements of the constant vector 
\fi\ << 1. Results in Fig. 4a show that indeed with this special choice of 
the perturbation, Axi(t) cx /i < Xi{t)xi{0) > +J2 < a;i(t)a;2(0) >, (noticing 
< a;j >= 0), in the linear response regime. This result is striking given that the 
perturbation SG is highly nonlinear. The response function shows damped os- 
cillations. Therefore, one expects Van Vleck-Weisskopf-Frohlich type resonance 
absorption [24]. The resonance is confirmed in Fig. 4b. Here we talk about the 
mapping Hamiltonian. The mapping "energy" is not the same as the physical 
energy. Nevertheless, one expects that the resonance absorption may manifest 
itself in the dynamics of x under periodic perturbation. We added a pertur- 
bation energy term 0.005 cos(a;i)(a;i -f- 2:2). Our stochastic simulations did find 
that the amplitude of the < Xi > oscillation driven by the perturbation shows 
a maximum near that suggested by the absorption spectrum (see Fig. 4c). 

3 Discussions and concluding remarks 

In this work, we showed how a recently discovered mapping between dissipative 
and Hamiltonian systems can lead to a general F-D relations. Further studies 
are needed for potential applications of the generalized F-D relation on study- 
ing systems described by Eqn [1] [12l [30] . As an example, one can use it on 
biological network reconstruction. The matrix M contains information on the 
network topology and parameter values. Researchers proposed to use the relax- 
ation dynamics to derive information of a biological network, especially a linear 
network |311 132| . A potential problem is that the relaxation dynamics alone may 
contain insufficient information to resolve a network topology and parameters, 
and to distinguish competing models. The F-D relation provides additional con- 
straints, and may serve as criteria for evaluating different models. It remains 
to be examined if one can distinguish abnormal (e.g. cancer) cells from normal 
cells based on their differences of relaxation and fluctuation dynamics. 



6 



We suggest that one should further exploit the analogy between the dissipa- 
tive and Hamiltonian systems. Many existing results for the latter may lead to 
new understanding for the former, provided appropriate quantities can be iden- 
tified. As an example, researchers have used oscillating signals to perturb bio- 
logical systems. They found that the oscillation amplitude of the variable being 
driven shows a maximum at certain driving frequency [331 134) . Here we showed 
that this phenomenon can be interrpreted as the Van Vleck-Weisskopf-Frohlich 
type resonance absorption [23]. The analogy may suggest further technical 
development in biological network studies parallel to the linear and nonlinear 
spectral methods used in physics and chemistry. Study analogous to Onsager's 
reciprocal relations is another example [51 |B] . A potential difficulty is to relate 
the mapped quantities to physically measurable quantities. 

The relaxation dynamics is important for characterizing the dynamics of 
a system. For example, the resistor model suggests that a HIV slows down 
its relaxation from an unstable active excited state to a stable dormant state 
through a series of intermediate steps in its regulation network [35) . Wang and 
coworkers found that many biological networks have rugged potential profiles 
(as defined by Eqn. [5]), which resemble what observed in protein folding stud- 
ies [3ni ISZl ISE] ■ The finding reinforces that dynamics of biological networks is 
complex. Implications of noises in biological systems have been extensively dis- 
cussed [33] ■ Here the mapping relates noise intensity to an effective temperature. 
Consequently we suggest that noise prevents a system from trapping into large 
numbers of undesired intermediate steady-states a complex system may have. 
Therefore, counter-intuitively existence of proper amount of noises, together 
with funneled landscapes [S^ , enhances rather than destroy the robustness of a 
network. Glass transition is an active research topic for the relaxation process 
of a physical system to an equilibrium state [301 Ell HH ■ A typical biological 
network is highly inhomogeneous, and is full of competing interactions. Thus 
its dynamics in some sense resembles that of spin glasses but with nonrandom 
quenching [35]. It is interesting to notice that the mapping between a dissipa- 
tive and a Hamiltonian system suggest possible glassy behaviors for relaxation 
to a nonequilibrium steady-state. If being identified, (noticing glassy state has 
been identified in protein folding studies [32 133]) also refer chapter 6 of [32] 
for a spin glass model in the form of Eqn. [Ij, it may have profound practical 
implications. For example, a patient might then have prolonged recovery period 
to the normal homeostatic state after a disease. One solution is to increase the 
temperature {i.e., noise intensity in the related regulation network) above the 
critical glass transition temperature. Indeed cells can regulate their noise levels 

In this work we only established a framework. Most of our discussions focus 
on biological examples. They equally apply to problems in other fields sharing 
the same mathematical structure, e.g., a stock market model described by Eqn. 

m 
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4 Appendix 



4.1 A 



Here we give the analytical solutions of the first example in case 1. The formal 
solution of the stochastic differential equations is, 



x(t) = < X > +exp(Kt)(x(0 

+z [ exp(K(i - T))C(T)dT 

-'0 



< X >) 



(19) 



where < x > is the solution of K • x + b = 0. The spontaneous fluctuation 
correlation functions are given by (with (3 — 1), 

{C^J{t)) {X^{t)- < X, >){Xj{0)- < Xj >) > 

cx j dx(0) cxp[-(x(0)- < X >)^ • • K • (x(0) 
- < x >)] (exp(Kt)(x(0)- < X >))^ {xj{Q)- < x, >) 



MK 



kj 



(1 + 4,) 



(20) 



The relaxation functions for the system relaxing from a steady state correspond- 
ing to fco + 6k to that of /cq ^^re given by, 



Axi{t) — Xi{t)— < Xi >cx 

J dx(0)exp[-(x(0)- < x' >)'^ 

(x(0)- < x' >)] (exp(Kt)(x(0)- 
cx (exp {Kt){< x' > - < X >)) . 



-M K 

2 

<x>)), 



0. 



(21) 



where < x' > satisfies K- < x' > + (fco + Sk, 0, 0) 

We calculate the transformation matrix M using the result of Kwon et al. 
[25] ■ For the choice of parameter values in this work, expression of the trans- 
formation matrix is given by. 



M = 



81z + 50z3 



81 + 25z2 -36 + 15Z -27- 20z 
36+l5z 81 + 9z2 45-12Z 
27- 20Z -45-12Z 81-1- IGz^ 



we want to point out that it is well known that the Langevin-type equations 
are only mathematical approximations of the physical processes, which may 
lead to unphysical results. For example, the concentration of a specie in a 
network may go to negative values. In the above derivations we assume that the 
probability of obtaining zero and further negative concentrations is negligible, 
and thus we can perform the integration from —oo to oo. 
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4.2 B 



For example 2 in case 1, the Langevin equations are propagated by 

x^itN) = x^{tN-i) + AtG,(x) + v/2zAV/3G(i), (22) 

where C,i{t) is generated from a Gaussian distribution with zero mean and unit 
variance, and we have set /3 = 1 throughout the work. We used At = 0.005 
in ah calculations. To calculate the relaxation function, we first propagated 
the network with a perturbation ( + 5k for results in Fig. 2, and vi + 5vi 
in Fig. 3) 4 X 10^ steps, so the system reaches a steady-state corresponding 
to the perturbed feg. Then we switch back to the u.nperturbed system, set 
time 0, record the values of Xi as Xi{Q), propagate for another 4 x 10^ steps. 
We calculated Xi{t) by averaging over 5 x 10^ such trajectories. The system 
is expanded linearly around the steady-state of the corresponding deterministic 
model. The predicted relaxation functions in Fig. 2 and 3 are calculated with 
5k — 0.1, and 5vi =0.1, respectively. We found that the system dynamics is 
within linear response regime when the predicted relaxation dynamics agrees 
with the simulated ones. 

For the example in case 4, the Langevin equations are propagated similarly. 
Results are averaged over 5 x 10^ trajectories. Rigorously speaking, there is a 
correction term to the stochastic Langevin equation. Our numerical algorithm 
adopts the Itoh interpretation. In [17], we discussed the relation between the 
Itoh interpretation and the zero-mass interpretation adopted for deriving the 
mapping Hamiltonian. The expression of the correction term can be derived 
straightforwardly. However, the correction term is proportional to 1//3, and is 
small for this system. Therefore we neglected this term here. 



4.3 C 

OA 



Here we prove that there is no need to consider Let's express pss — ■N'po, 
and (l){x) = — 4(ln(0o + lnA/'+ (t>o), then. 



_,[<X«,M)>-<X><||>1 

=< xio "'"";'"'"'' ) > - < X X ^ >i 

OA OA 



We thank Kenneth S Kim for discussions. 
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Figure 1: Numerical tests with a linear network, (a) The network being tested, 
fco = fci = fc2 = fcs = ^4 = 1. (b) The normalized relaxation function Axi{t) = 
{xi{t)- < xi >) /(xi(0)- < xi >). (c) The Ax2{t) calculated with Eqn [TTl 
overlap with the exact result from Eqn [^l] (solid line) , but that calculated from 
Eqn [12] with M = 1 gives wrong prediction (dashed line). Similar results for 
Ax^it). 
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Figure 2: Reaction network under Michaelis-Menten kinetics. Model parameters 

ko = l,Vi = l,Ki = 0.5 {i = 1 — 4), except V2 = 1.5, z = 0.005. (a): The steady- 
state distribution for the unperturbed system, (b)-(d): The relaxation function 
with various perturbation 6k. To compare results with different 6k, all results 
are normalized by the maximum of |Aa;t(t)|. 
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Figure 3: Same as Fig. 2, except the perturbation is on vi. 
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Figure 4: Numerical tests with a model showing limiting cyele dynamics, (a) 
The normalized relaxation function Axi{t) = {xx{t)— < Xi >) /(xi(0)— < Xi > 
) from simulation and prediction of the F-D relation, fi = = 0.005. (b) 
The "absorption spectrum" calculated from the fluctuation correlation func- 
tions with fi = f2 = 0.005 cos(wt) (c) The simulated amplitude of the < Xi > 
oscillation driven by an oscillating perturbation used in b. 
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